function X_peaks = ecg_peaks(X,xlimlow,xlimup,Thres1,Thres2,type)
% ECG_PEAKS         ECG Peak detection
% X_peaks = ecg_peaks(X)
% X_peaks = ecg_peaks(X, Thres1)
% X_peaks = ecg_peaks(X, Thres1, Thres2, type)
%
% This function detects position of the QRS complex. First, signal is 
% smoothed by ecg_smooth(). Next, turning points are found by 1st order 
% difference analysis. Then 'smooth' turning points are removed by 2nd 
% order difference analysis, leaving only 'sharp' peaks. 
%
% Inputs:
% X: discrete time signal
% Thres1: Relative threshold in 1st difference analysis (0~1). 
%         The higher the value, the more data points are selected.
% Thres2: Relative threshold in 2nd difference analysis (0~1). 
%         The higher the value, the less data points are selected.
% Type: 'impulse' or 'window'
%
%
% Copyright: Group 8
% Modified at 20th Oct 2010
%-------------------------------------------------------------------------
if(nargin < 6)
    type='impulse';
end
if(nargin < 5)
    Thres2=2.3e-3;
end
if(nargin < 4)
    Thres1=0.01;
end
if nargin < 3
    xlimup = 13;
    xlimlow = 10;
end
%Any number is smaller than eps is considered to be zero
eps1 = max(abs(X))*Thres1;
eps2 = max(abs(X))*Thres2;

%Smooth the signal without shifting the peaks
X_smooth = ecg_smooth(X,9,101);

%First pass: find all diff(X,1)==0
X_pass1 = find_zeros(X_smooth,1,eps1);

%Second pass: find all abs(diff(X,2))<eps2
X_pass2 = find_zeros(X_smooth,2,eps2);

%Find only 'sharp' peaks
X_windows = (X_pass1>X_pass2)';

if(strcmp(type,'impulse'))
%Change windows to impulses
X_peaks(2:length(X_windows)) = diff(X_windows)>0;
X_peaks(1)=0;
elseif(strcmp(type,'window'))
X_peaks = X_windows;
else
error('Wrong type: only "window" or "impulse" are suppported')
end

%Debug mode: draw graph if no output is specified
if(~nargout)

t=((0:length(X)-1)/550)';

%Window obtained form first difference analysis
t_p1 = t(X_pass1);
X_p1 = X(X_pass1);

%Window obtained form second difference analysis
t_p2 = t(X_pass2);
X_p2 = X(X_pass2);

%Final selected impulses / windows
t_selected = t(X_peaks);
X_selected = X(X_peaks);

plot(t,X,'r',...
	t_selected,X_selected,'xg','markersize',15,'linewidth',3);
legend('signal','peaks');
xlim([xlimlow xlimup]);
end


function X_zeros = find_zeros(X,N,eps1)
if(N<=0)
    N=1;
end
len = length(X);
diff_X = diff(X,N);
X_zeros((N+1):len) = abs(diff_X) < eps1;


